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Abstract. 

We study the impact of neutrino masses on the shape and height of the BAO peak 
of the matter correlation function, both in real and redshift space. In order to describe 
the nonlinear evolution of the BAO peak we run N-body simulations and compare 
them with simple analytic formulae. We show that the evolution with redshift of the 
correlation function and its dependence on the neutrino masses is well reproduced 
in a simplified version of the Zel’dovich approximation, in which the mode-coupling 
contribution to the power spectrum is neglected. While in linear theory the BAO peak 
decreases for increasing neutrino masses, the effect of nonlinear structure formation 
goes in the opposite direction, since the peak broadening by large scale flows is less 
effective. As a result of this combined effect, the peak decreases by ~ 0.6 % for 
= 0.15 eV and increases by ~ 1.2% for = 0.3 eV, with respect to a 

massless neutrino cosmology with equal value of the other cosmological parameters. We 
extend our analysis to redshift space and to halos, and confirm the agreement between 
simulations and the analytic formulae. We argue that all analytical approaches having 
the Zel’dovich propagator in their lowest order approximation should give comparable 
performances, irrespectively to their formulation in Lagrangian or in Eulerian space. 
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The baryon acoustic oscillations (BAO) peak of the galaxy correlation function 
(CF) is one of the prominent observables in present day cosmology [H [2]. The peak 
position corresponds to the acoustic scale at decoupling and can be used to infer the 
angular diameter distance as a function of redshift and the Hubble parameter H{z). It is 
well known that nonlinear effects of structure formation cause a broadening of the peak, 
which degrades the measurement of the acoustic scale by roughly a factor of 3 at ^ = 0 
[3]. The blurring of the peak is due mainly to the random displacement of the galaxies 
from their original location, by an average of ~ 6Mpc/h. Reconstruction techniques 
have been developed H El El m El [9] to undo the effect of these bulk motions and get 
back a CF closer to the linear one, and have been successfully applied to real data (see 
for instance [IQ]). The effectiveness of these techniques, which are based essentially on 
the Zel’dovich approximation or on variants of it, can be seen as a proof, a posteriori, 
that the physics of the BAO peak degradation is well understood and under control. 

In this paper we explore the possibility of extracting cosmological information from 
these nonlinear effects, rather than seeing them as a noise to be removed. In practice, 
we will shift the focus from the peak position to the peak shape of the nonlinear CF, 
and we will study the latter as a function of redshift and of cosmological parameters, in 
particular neutrino masses. 

Measuring the total mass of neutrinos by using cosmological data set is one of the 
most important goals of present and future large scale structure surveys. At present 
the tightest constraints are obtained by adding to the cosmic microwave background 
data, BAOs and the Lyman—a forest and result in a 2cr C.L. upper limit of 0.14 eV 
mm- While BAOs are helpful in breaking degeneracies between the cosmological 
parameters probed by the cosmic microwave background (namely and Hq), the 
Lyman-a forest data or other small scale structure probes like weak lensing and galaxy 
cluster data probe the neutrino induced suppression on the matter power spectrum. 
Although non zero neutrino masses have been claimed in the recent past in order 
to reconcile putative tensions between cosmic microwave background data and lower 
redshift probes, it appears that this solution is not supported by BAOs and the recent 
Planck results [nmi]. In any case the tight upper limits above, if conhrmed, will have 
tremendous implications for particle physics experiments [T^ and future cosmological 
searches like Euclic^ These new Endings are thus providing a strong motivitation to 
further explore neutrino induced non-linearities in unexplored regimes, like the one that 
we will be addressing in this paper. 

Cosmologies with massive and massless neutrinos sharing the value of cr8,ch[§] will 
have very similar halo mass functions [H]. This implies that constraints on the sum 
of the neutrino masses from observables such as cluster number counts will be very 
degenerate with ag (the well known degeneracy) |T7]. One way to break that 

degeneracy is to use observables that prove larger scales, like the anisotropies in the 

f http://www.euclid-ec.org 

§ This is the value of (Ts computed using the CDM plus baryons power spectrum; in the following, 
tim = 1 — Pa denotes the present total dark matter density. 
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Figure 1. Nonlinear matter CFs in real space and at redshift z = 0. The solid curves 
are for the massive neutrino cosmologies. The dashed curves are for cosmologies with 
massless neutrinos, but with the same values of tlb, fim, and crs.cb as the corresponding 
cosmology with massive neutrinos. The CF has been computed as described in the 
text, and here and in all the other plots in this paper, it has been multiplied by i?^. 


cosmic microwave background or the galaxy correlation function on the BAO scale. In 
Figure we show the nonlinear CF for massive neutrino cosmologies (solid lines) and 
corresponding massless neutrino ones with the same cr8,cb, k2h, and Qm (dashed lines). 
As we see, the BAO peak provides a clearly distinct signature for these cosmologies. 

We will show that a very simple approximation to the nonlinear power spectrum 
(PS), consisting in removing the mode-coupling part from the Zel’dovich approximation, 
gives an accurate description of the BAO peak in the CF, and an even better one of 
the ratio between CFs of different redshifts or of cosmologies with different neutrino 
masses. We also explore improved approximations, which give a better PS and CF, but 
do not improve substantially the CF ratios. Of course, the mode-coupling part is very 
important to recover the full nonlinear PS, where it becomes dominant at small scales, 
but since it is quite smooth on the BAO scales its contribution to the CF is limited to 
a subpercent shift in the acoustic scale, and almost disappears in CF ratios. 

We will consider the CF both for matter and halos, in real and in redshift space, 
and compare our analytic formulas to N-body simulations, finding perfect compatibility 
in all cases. 

The paper is organised as follows: in Section we discuss our analytic 
approximation to the CF in the BAO region, both in real and in redshift space, first for 
matter and then for halos; in Section we present our suite of N-body simulations for 
ACDM and for massive neutrino cosmologies, and describe the algorithm to measure 
the CF; in Section we compare our analytic formulas with the CF measured from our 
simulations and also, for the ACDM cosmology, with the FrankenEmu [TS] N-body based 
emulator; in Section we summarize our results; in Appendix A we derive the general 
expression for the nonlinear PS, eq. ([^; in Appendix B we show the relation between our 
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simplest approximation of the nonlinear CF, eq. Q, and the Zel’dovich approximation; 


m 


Appendix C| we derive the PS in redshift space inclnding the effect of bnlk flows. 


Finally, in Appendix D we discuss the dependence of our N-body simulations on mass 
and force resolution. 


1. Nonlinear evolution of the BAO peak 

In full generality, the nonlinear matter PS at redshift z has the following structure 

[iniED]: 

P(fc, Z) = G\k, z)P^^^{k, Zin) + PMc{k, Z) , (1) 

where P''^'^{k,z) = D‘^{z)P^'^"‘{k, Zin) is the linear PS, and Zin is some initial redshift 
chosen well after decoupling and such that all the relevant scales are still in the linear 
regime (it can coincide with the redshift at which we start the N-body simulations, 
which in this paper is Zin = 99, see Section]^. The nonlinear effects are completely 
encoded in the two functions appearing at the RHS: the propagator G{k, z), representing 
the cross-correlator between the nonlinear density field at redshift z and the initial one 
at Zin [IS El], and the “mode-coupling” term PMc(,k,z). 

We stress that the above expression is completely general, the only assumption 
behind it being that the nonlinear density field, (5(k, z), is some “functional” of the 
initial density and velocity fields, see [Appendix A| for details. Then, one can compute 
these quantities in any consistent approximation scheme, such as Eulerian or Lagrangian 
perturbation theory (PT). 

The BAO wiggles of P^^^{k,z) are in general smoothed out in PMc{k,z) as its 
computation involves momentum integrals in which two or more linear PS evaluated 
at different scales are convolved. Therefore, as we will demonstrate below, the BAO 
information is basically confined to the G'^{k] z)P^^^{k) term which, after Fourier 
transform, accounts for approximately all the BAO peak in the CF. 

The propagator has been studied thoroughly in the recent literature 
In Zel’dovich approximation it is given by 


see 


Appendix B) 


G^^^%k,z) =e-F 


( 2 ) 


where (t‘^{z) is the 1-dimensional velocity dispersion evaluated in linear theory, namely. 




Pq P^^‘^{q,z) 


(27r)3 


( 3 ) 


In the following, we will mostly consider the CF obtained in Zel’dovich approximation 


by neglecting the mode-coupling part (see eq. (B.7)). Namely, we will study the Fourier 
transform of 


P^^\k,z) = 


( 4 ) 


and we will show that it gives a very good approximation to the nonlinear matter CF 
in the BAO peak region, and in particular to ratios of CF’s for different redshifts, or 
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for different values of neutrino masses, both in real and in redshift space. We will also 
show how it can be extended to describe halos. 

In Figure 1^ we plot the linear CF, the CF obtained from Q, and the CF from 
our N-body simulations for a ACDM cosmology (the corresponding list of cosmological 
parameters is given in Section]^. We also plot the difference between the N-body based 
interpolator FrankenEmu [18] and the CF obtained from (|^, namely the CF associated 
to 

P{k,z)-P^^\k,z) = PMc{k,z) + (5) 

As we anticipated, the simple formula (|^ captures most of the physics responsible 
for the BAO peak broadening (he. the difference between linear and nonlinear BAO 
peaks), while neglecting the fully nonlinear mode-coupling term and approximating the 
full propagator with its Zel’dovich limit gives negligible and mostly scale-independent 
contributions at these scales. 

As it has been emphasised in [211 |2S|, coherent bulk flows on scales larger than the 
BAO scale Ibao — 110 Mpc/h cannot have any effect on the BAO peak in the CF. This 
is a manifestation of the general property of Extended Galilean Invariance (EGI) or of 
the Equivalence Principle, discussed in [26l EZl EHl EHl EQ]. Therefore, one would expect 
that the effect of modes q ~ 271/Ibao should be cut-off in the momentum integral in 
(|^, in order to compute the GF at scales 0{Ibao) and smaller. On the other hand, 
the full expression ([^ automatically exhibits EGI: the decoupling of the large modes is 
achieved via a cancellation of their effect between the propagator term, without any IR 
cutoff in the momentum integral, and the mode-coupling one. The point is that, as we 
have seen, the mode-coupling term is smooth on the BAO peak scale, and therefore one 
can, as a first approximation, neglect it if one is interested in this feature. 

This approximation can be improved in two respects. On the one hand, the 
Zel’dovich approximation uses the linear expression for the displacement field at all 
scales, and will therefore fail for wavenumbers q larger than some value A in ([^. 
On the other hand, eq. (|^ contains no mode-coupling effect at all, while these are 
known to induce some shift in the BAO peak position, although at a subpercent level 
[3ll[l9l|32l[33l|3l|. 

We can consider an improved version of our formula by using standard FT (|^, 
namely 

P^‘^\k,z) = e-'’'<Ab) (p^-(fc,z) + P22 (A:,z) + 2P'*’^(fc,z)EA,oo(A;,^)) , (6) 

where the 1-loop mode-coupling term is given by 


P22{k,z) = 




98 (27r) 


dqP^^^{q,z) / dxP^^^{^/k^ + q^ -2qkx,z) 




'-1 






( 7 ) 
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while 

roo 


k + q 
k — q 

( 8 ) 

and o'^^jy(z) is given by eq. ^ in which a UV cntoff on the momentum q at the scale A 
has been imposed. Using the fact that for soft modes k^a^j^{z) ~ —2I]Q^\{k,z), we see 
that expanding the exponential in (|^ one recovers the 1-loop expression for the PS in 
standard PT, 

P^^^k, z) = P^^^ik, z) + PisiK z) + P 22 {k, z) + 0(2 loop). (9) 

The extension to higher loop orders, possibly including also a better treatment of 
the short modes along the effective approaches discussed in |35l IMl 125], is possible, 
although computationally more and more demanding. 

Our main point in this paper is to show that the physics of the BAO peak 
degradation is well understood, so that by measuring not only the peak position, but 
also the peak shape, one can recover cosmological information. To discuss this point 
we will focus on the simplest approximation, eq. (Q, keeping in mind that it can be 
systematically improved. 


158 + 100 ^ 




q'^ ^ 3k'^ I'q^ 

q 


3 VP 


2 + 7^ ) log 



1.1. Parametric dependence of the BAO peak 


We will concentrate on the BAO peak in the CF. We start by discussing pure matter in 
real space, for which the CF is given by 

1 

^(72) = y dqqAn{qR)P{q) , (10) 

where we have omitted the time (or redshift) dependence. We will also use the 
“moments” ^ 1 ^ 2 ( 7 ?), defined as 

in{R) = dqq {qRY sin(gi?) P{q). (11) 

In the BAO peak, R = R, defined by 


e'(^) = 


d^R) 

dR 


= 0 , 


( 12 ) 


we have ^ dqq^ cos {qR) P {q) = f, (7?), and so ^"(7?) = —^2{R)/R^- 

If we take the simplest expression (|^ for P{q), then the dependence of the CF ^{R) 


on af is given by 


dlog^jR) 

dlogaf; 


7?2 e(i?) 


( 13 ) 
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If we consider two linear PS differing by and, correspondingly, two af 

differing by 


3 

the CF changes by 

sm = 


dfq 6P’'^‘^{q) 
(27r)3 


(14) 


27r^R 


dqq sm{qR) 6P’'^"'{q) e 


6 {R) da] 




1 


= ^ I dqq" dP^^'^iq) 


/ sm{qR) 2^2 _ 1 6(^) \ 
\ qR 3 qm^ ) 


lin 


(15) 
(9) = 


(16) 


For instance, if the two linear PS differ only in their normalizations, hP 
P’‘^‘^{q)6A/A, the corresponding height of the BAO peak changes by 

dAm ^d_Af alUR) \ 

m A ^ R^m J ’ 

where the second term inside parentheses gives the -scale-dependent- nonlinear effect. 
A special case is that of considering two different (and nearby) redshifts, in which case 
5A/A = 6D{zy/D{z)‘^. In this case, at z = 0, the nonlinear term gives typically 
a 15 — 25% negative correction at the peak position and a 20% positive one at the 
minimnm left to the BAO peak, with respect to the linear resnlt. 

Eq. (15) can also be nsed to estimate the peak change between two slightly different 
cosmologies. For instance, let ns consider the change of the BAO peak at z = 0 between 
two cosmologies with different nentrino masses (the two cosmologies only differ from each 
other by the nentrino mass, and by the cold dark matter abnndance, in snch a way that 
the total matter component, Om, is the same). For '^m^, = 0.15 eV, eq. (15) predicts 
a decrease of the peak height of ~ —0.6% with respect to the massless nentrino case. 
At the higher masses that we have considered the peak heights is instead greater than 
in the massless case. Specihcally, eq. (15) gives an increase of ~ 1.2% for = 0.3 

eV, and of ~ 5.7% in the case of = 0.6 eV. This behavior is conhrmed by the 

resnlts presented in Fignrej^ In particniar, we see that nonlinear effects invert the trend 
with respect to linear theory. Indeed, in linear theory the BAO peak height decreases 
with increasing nentrino masses (see the red-dashed cnrves at i? ~ lOOMpc/h in the 
hgnre). The bnlk flows, on the other hand, are less effective for higher nentrino masses 
(as the corresponding af^ is lower) in degrading the linear BAO peak, and therefore 
they increase the ratio of the height in the massive vs. massless case with respect to 
the linear prediction. For Jfmiy = 0.15 eV, the decrease of the peak predicted by the 
linear theory dominates over the increase dne to the bnlk flows. The opposite is trne 
for = 0.3, 0.6 eV. 

In Section we compnte the CF in real and redshift space, starting from the power 
spectra introdnced above. Specihcally, in real space we compnte and show 

1 

dq q sin {qR) (g) , 


e^{R) = 


2ti^R 
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1 

‘'■(-)[P””(<,) + F22(<;)] . (17) 

where the last one is obtained by taking the A —)■ cx) limit in eq. §. As we will see, 
reproduces the shape and position of the BAO peak slightly better than However, 
they perform equally well in reproducing ratios of the CF between different redshifts 
(within a fix cosmology) or ratios of correlations functions of different cosmologies (at 
equal redshift). For this reason, we only use in the subsequent analysis. 

Specifically, after studying the matter CF in real space, we study the matter CF in 
redshift space, and the CF for biased halos that are predicted from and we compare 
them with the results of the numerical simulations. We verify that the agreement is 
comparable with that obtained for the matter CF in real space. 


1.2. Redshift space 


In redshift space, the PS is a function of time, of the magnitude k of the momentum of 
the modes, and of the cosine /r = k ■ z of the angle between the line of sight (that we 
take along the z-axis) and k. On large scales the PS is given by the Kaiser formula 


(k, IJ.) ^ (1 + l/ fY p"’- (k) , (18) 

where / = and a is the scale factor (to simplify the notation, the time dependence 

of / is understood; in the following, from now on, also the time dependence of cx^ is 
understood). The same effect on the nonlinear PS (|^ gives 

pW (A), fj) = (A), /i) ^ 

Appendix C 


see 


The CF has also nontrivial angular dependence, fg (R). In the following, we only 
study its monopole component. For a generic PS that is function of k and pf, the 
angular averaged CF is 


6(R) = - / dnf,^g{R) = 


2ti^R 


dq q sin (qR) / dp Pg {k, p) 


( 20 ) 


Using the approximations (18) and (19) in (20) we get, respectively. 


/kaiser , 




(p)= 1 + -/+V r(^), 


{R) = 


2'k‘^R 


dqqsiYi{qR)P^^{q)e-‘^"^"p\^2V]-qa,, /, 1 


( 21 ) 


where we have defined 


t 3 + 417^71 + 7100 0FErf[/f/7] 3+2(2A + 

’P )- « ■( ) 
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We are also interested in CFs of dark matter halos. Following [371 |38l ES] , we consider 
the following large scale bias predicted by the peaks model between the halo and the 
matter density helds 

6h{k,z)) 


= [bio {z) + boi {z) k"^] 




(23) 


{b'ff{k,z)6^ff{k,z)) 

where, in addition to what done in [3Zl ESI EH], we have added the exponential 
factor to account for the (unbiased) bulk flow suppression. We extract the coefficients 
^105 ^01 by matching at largest scales the correlation function obtained from our N-body 
simulations. 

From this relation, we obtain the real space halo-halo CF 


hi’ (fl) = 


2tt^R 


dq q sin {qR) [6io -h &oi (f] e 


(24) 


Strictly speaking, this relation is only valid on large scales, but we show that 
assuming it to be valid at all scales allows for a good reconstruction of the BAO peak 
in the halo CF: the measurement of the halo CF from the N body simulations has a 
greater uncertainty than that of matter, due to the smaller statistics (there are fewer 
halos than dark matter particles in the simulation; this is particularly true at higher 
redshift, where the halo data become progressively more uncertain). We hnd that the 
CF (24) agrees with the numerical simulations as well as in the matter case, compatible 
with the increased error bars. 

We notice that using N-body simulations it was found by [iQl |16] that the halo 
bias in cosmologies with massive neutrinos is scale dependent, even on very large scales. 
The reason for this is that dark matter halos are biased tracers on the underlying 
CDM-|-baryons held, not of the total matter one. If the halo bias is dehned with respect 
to the underlying CDM-|-baryon held, the halo bias turns out to be scale independent 
on large scales and universal, as found by [16]. Therefore, in cosmologies with massive 
neutrinos it is better to compute the halos correlation function using the above equation 
but employing the properties of the CDM-|-baryons held rather than these of the total 
matter. Given the fact that the error bars from our halo catalogue are very large, 
and since the scale-dependence bias ehect induced by massive neutrinos is small, we 
decided, for simplicity and clearness, to employ the above equation even for cosmologies 
with massive neutrinos. 


Finally, we are interested in the halo CF in redshift space. We start from (24), and 
assume that the halo velocity is unbiased, so to obtain the redshift space PS 


Using (20), we obtain 

^ilh (R) = 


dq q sin (qR) P (q) e 


lin „-q^cr'^ 


R 


\/2 + / qcy, f,bio + q'^boi 


( 26 ) 


where R has been dehned in (22). 
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We have run N-body simulations for four different cosmological models; one with 
massless neutrinos and three with massive neutrinos: = 0.15 eV, = 0.3 

eV and ~ O-® simulations share the value of the following cosmological 

parameters, which are in excellent agreement with the latest results by the Planck 
collaboration [13]: flm = 0.3175, fib = 0.049, IIa = 0.6825, h = 0.6711, Ug = 0.9624, 
As = 2.13 X 10“®. In cosmologies with massive neutrino we set the value of flcdm to 
flm — fib — fli/, where = ^mjy/(94.1 eV). The value of ug obtained from these 

parameters is 0.834, 0.801, 0.764, 0.693 for the models with Y'^f = 0-0; 0-3 and 

0.6 eV neutrinos, respectively. 

For each cosmology we have run 100 independent N-body simulations [f] using the 
TreePM code GADGET-III [H]. In each realization we follow the evolution of 256^ 
cold dark matter (plus 256^ neutrino particles for massive neutrino models) down to 
z = 0 in a periodic box of 1000 Mpc/h length. The Plummer equivalent gravitational 
softening is set, for each particle type, to 1/40 of the mean inter-particle linear spacing. 
We save snapshots at redshifts 2, 1, 0.5 and 0. 

We have generated the N-body initial conditions, at Zin = 99, in the following way: 
we set the particles in a regular cubic grid and then we displace their positions according 
to the Zel’dovich approximation. The linear matter power spectra and transfer functions 
are computed using the GAME code [l^. In simulations with massive neutrinos we use 
a transfer function which is a mass-weighted average between the CDM and baryons 
transfer functions to set up the initial positions and velocities of the CDM particles. 

Dark matter halos are identihed using the Friends-of-Friends algorithm [13] taking 
a value for the linking length parameter equal to b = 0.2. We only consider halos 
containing at least 32 CDM particles, i.e. halos with masses larger than ~ 1.6 x 
10^^ h“^M 0 The two-point correlation function of the dark matter halos is computed 
using the Landy-Szalay estimator H 

DD{r) — 2DR{r) + RR{r) 


^(r) = 


RR{r) 


(27) 


with DD{r) and RR{r) being the number of pairs with distances in the interval 
[r — Ar, r -|- Ar] in the halo and random catalogues, respectively. DR{r) is the number 
of halo-random points pairs in the above interval. Our random catalogue contains half 
million points, which is a number more than 30 times larger than the typical number of 
halos in one realization at x = 0. 

Unfortunately, the use of the Landy-Szalay estimator to compute the two-point 
correlation function of the matter held is computationally unfeasible, given the large 
number of tracers in our simulations. Therefore, we evaluate the matter correlation 
function using the estimator presented in [15] that we now briehy describe. We hrst 
evaluate the density contrast held, S{r), using a grid with 256^ (or 512^) cells using the 

t We associate an error to the mean two-point correlation function equal to the variance around the 
mean from the 100 realizations run for every cosmological model divided by vTO. 
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cloud-in-cell (CIC) mass assignment scheme. We then Fourier transform (5(r) to obtain 
5(k) and compute in each point of the grid the quantity |(5(k)p. Finally, we Fourier 
transform the previous quantity and compute the two-point correlation function by 
averaging over all modes that he within the same r—bin. We set the width of the bins 
in radius equal to where L is the size of the simulation box. 

In order to verify the robustness of the above algorithm we have also computed the 
two-point correlation function using the estimator presented in [IHlISIIIH]- [|] We find an 
excellent agreement among the correlation functions computed using the two different 
estimators. We notice that matter correlation function results shown in the paper have 
been obtained using the former estimator. 

Correlation functions in redshift-space are computed using the above methods after 
displacing the particle (or dark matter halo centers) positions from real-space to redshift- 
space using eq. (C.2). We have created three different catalogues, by displacing the 
particles/halos along each cartesian axis according to the peculiar velocities projected 
along that axis, and in the paper we show the results when averaging among the three. 


3. Results 

In this section we discuss the dependence of the CF on redshift and on values of the 
neutrino masses, and compare the results of N-body simulations to the analytic formulas 
discussed in Sect. [H 

In Figure]^ we show the matter CF in the massless neutrino case at z = 0 (first 
row) and at z = 1 (second row). The right panels are a zoom of the left panels centered 
at the BAO peak. We see that linear theory predicts a substantially higher BAO peak, 
and we clearly see the nonlinear broadening in and In the left panels of the 
figure, we also show with a black solid (dashed) line the residual difference between the 
CF obtained from the FrankenEmu [18] N-body based emulator, and (and 
This difference - which for is mostly due to the mode-coupling contribution to the 
PS - appears to be very weakly dependent on R and on the redshift. We stress that the 
exponential factor in Q is not fitted to the data, but it is obtained from the linear PS, 
according to ([^. 

We also note that the addition of the 1—loop P 22 term in shifts the peak to 
slightly smaller scales, at the subpercent level pH[T9l|32l|33l|3l| , and it slightly improves 
the agreement with the FrankenEmu result at z = 0 (the difference plotted as a black 
dashed line is both slightly closer to zero, and slightly more scale independent than the 
difference plotted as a black solid line; this is much less marked at the higher redshift 
shown, since the 1—loop P 22 term becomes more irrelevant at increasing redshift). 
However, as we show in the next figures, the addition of the 1—loop P 22 term does 

f We first compute the density field in a grid with 160^ cells using the nearest grid point 
mass assignment scheme and then we evaluate the correlation function using the estimator ^(r) = 
jA — W- ■ SiSj , where the sum is performed over cells whose distance lie in the interval \r — Ar, r-l-Arl 
and fVpairs is the number of cell pairs in that interval. 
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Figure 2. First row: Matter CF in real space, for massless neutrinos, and at redshift 
z = 0. The right panel is a zoom of the left panel centered at the BAO peak. The 
data points are from our N-body simulations; the red dashed, green solid, and blue 
solid lines are, respectively, ^^^\and^(^\ defined in eq. (17), multiplied by 
The black solid (dashed) line at small values in the left panel is the difference ([^ 
between the CF from the FrankenEmu m N-body based emulator and (and 
also rescaled by R^. The black solid line in the right panel is the FrankenEmu CF, 
times R^. Second row: same as in the first row, but at redshift z = 1. 



Figure 3. Ratio between the matter real space CF at two different redshifts, for 
massless neutrinos. The top, middle, and bottom curves in the figure are ratios of CF 
at z = 0.5, z = 1, z = 2, respectively, divided by the corresponding CF at z = 0. The 
data are ratios between our N-body simulations; the red dashed, green solid, and blue 
solid lines are ratios between, respectively, ^(^\and^^^\ defined in eq. (17). The 
black solid lines are ratios between CF obtained from the FrankenEmu emulator. 
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R [Mpc/h] R [Mpc/h] 


Figure 4. Matter CFs in real space. Analogous of Figures]^ {z = 0) and|^ (ratios 
between CFs at different z), but now for massive neutrinos. The figures in the first 
row are for ffrrii, = 0.15 eV, while those in the second row are for rrii, = 0.3 eV. 


not improve significantly the ratios between CF’s. 

Indeed, while the CF does not perfectly reproduce the N-body CF, it tracks 
extremely well how the CF changes with redshift. We see this from Figure where we 
show ratios between matter CF (of the same cosmology) computed at different redshift. 
The ratios obtained from are in excellent agreement with the ratios obtained from 
our N-body data, as well as with the FrankenEmu. We also see that, as we just 
mentioned, the inclusion of the P 22 term does not provide a signihcant improvement 
on these ratios. 

Identical conclusions are obtained in the comparison between and our N-body 
data in the case of massive neutrinos. Notice that FrankenEmu does not provide data 
for these cosmologies. We show this in Figure]^ where we present the CF at 2 ; = 0, and 
the ratio between CFs at different redshift, in the case of = 0.15 eV (hrst row) 

and 0.3 eV (second row). In these cases, we computed the velocity dispersion using 
the linear PS for total matter in eq. ([^, that is for Sm = ^c<^c + as it is 

the source of the Poisson equation. 

It is natural to ask whether an equivalent agreement takes place also in redshift 
space. This is conhrmed by Figure where we show the comparison between the 
angular-averaged redshift space CF (21) and the one obtained from the N-body data. 
The linear correlations functions in real and redshift space are related to each other by 
the Kaiser relation ( [2l| ). Not surprisingly, this also overpredicts the BAO peak. On the 
contrary, the CF .^i^^shows an equal agreement with the N-body simulations as its real 
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R [Mpc/h] R [Mpc/h] 


Figure 5. Matter CFs in redshift space at z = 0 (first column) and ratios between 
CFs at different redshift (second column). The figures in the first, second, third row 
are for ffrui, = 0, 0.15, and 0.3 eV, respectively. The data are from our N-body 
simulations; the red dashed, and green solid lines are, respectively, for and for 

defined in eq. (211. 


space counterpart 

The real space CF and its redshift space counterpart are an optimal tool to 
study the dependence of the CF on the neutrino masses. This can be seen from Figure 
1^ where we show ratios between the CF of a cosmology with massive neutrinos divided 
by a cosmology with massless neutrinos (the two cosmologies only differ from each other 
by the neutrino mass, and by the cold dark matter abundance, in such a way that flm is 
the same for them). Also in this case, the ratios obtained from (left column plots) 
and (right column plots) are in excellent agreement with the ratios obtained from 
the N-body data. 

Actual measurements of the BAO peak involve biased objects. In Figurej^we study 
the agreement between the halo CF (24) and the halo CF obtained from our N-body 













































The effect of massive neutrinos on the BAO peak 


15 



Figure 6. Ratio between the z = 0 matter redshift space CFs of two cosmologies 
with different neutrinos masses. The first, second, and third row show the CF for 
= 0.15, 0.3, 0.6 eV, respectively, divided by the corresponding CF for massless 
neutrinos. The left column shows the ratios in real space. The data are ratios between 
our N-body simulations; the red dashed, green solid, and blue solid lines are ratios 
between, respectively, ^*-^\and^^^\ defined in eq. (17). The right column shows 
the ratios in redshift space. The red dashed, and green solid lines are, respectively, 
Raiser ^(1)^ defined in eq. (Eil. 


simulations, as described in the previous section. The comparison is less probing than in 
the matter case, due to the increased sample variance of the latter (there are fewer halos 
than dark matter particles in the simulations). This is particularly true at increasing 
redshifts, and for this reason we only show halo data a.t z = 0, 0.5. The two solid 
lines shown in the figure are obtained with either a constant density bias, b (k) = bio 
(green line) or a bias of the type b {k) = 610 + &01 (blue line), times the exponential 
suppression due to the bulk flows, see eq. (24). The bias coefficients are obtained by 
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2 niy = 0.0 eV; (real space) 


R [Mpc/h] 


R [Mpc/h] 




R [Mpc/h] 

^ niv = 0.15 eV; (real space) 



R [Mpc/h] 




R [Mpc/h] 


R [Mpc/h] 



R [Mpc/h] 



R [Mpc/h] 



R [Mpc/h] 


Figure 7. Halo CFs. The figures in the first, second, and third row are for 
ffnii, = 0, 0.15, and 0.3 eV, respectively. The left panels show the real space CF 
at z = 0; the middle panels show the ratio between z = 0.5 and z = 0 CF in real space; 
the right panels show the same ratio in redshift space. The data are from our N-body 
simulations. The green (respectively, blue) curves are obtained from the bias function 
bio (respectively, biQ + 6oi fc^)- 


fitting the N-body correlation function (23) at large scales. For the plots on the left 
column we see that already using the constant bias allows to reproduce the height of 
the BAO peak. Allowing for the 6 oi term improves the agreement with the N-body 
CF at values of R smaller than the peak. The plots in the second and third column of 
the hgure show ratios of CF at different redshifts (the second column shows ratios in 
real space, while the third column shows ratios in redshift space). We see that both the 
ratios obtained from a constant or from a linear bias are in agreement with the ratios 
from the N-body data (in fact, not appreciable scale dependence can be observed from 
the N-body ratios due to their error bars). 

Finally, in Figure we show ratios between the halo CF at 2 ; = 0 of a cosmology 
with massive neutrinos and of a cosmology with massless neutrinos (keeping the same 
klm for all cosmologies, as we did for Figure [^. We see that the ratios obtained using 
the 610 + &01 bias are in better agreement with the central value of the N-body data 
than those obtained from a constant bias, although also the latter appear to be within 
the 1 — a error bars of the N-body data. 
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Figure 8. Ratio between the z = 0 halo redshift space CFs of two cosmologies 
with different neutrinos masses. The first, second, and third row show the CF for 
ffrifii, = 0.15, 0.3, 0.6 eV, respectively, divided by the corresponding CF for massless 
neutrinos. The left (respectively, right) column shows the ratios in real (respectively, 
redshift) space. The data are from our N-body simulations. The green (respectively, 
blue) curves are obtained from the bias function bio (respectively, 6io + 6oi)- 


4. Conclusions 

The shape of the BAO peak contains cosmological information. In this paper we have 
focnsed on nentrino masses, and we have shown that, comparing cosmologies with the 
same Qm, massive nentrinos with ruy = 0.15 eV decrease the peak height by ~ —0.6% 
with respect to the massless nentrino case, whereas = 0.3 eV (0.6 eV) increase 

it by 1.2% ( 5.7%). These effects shonld be measnrable in futnre snrveys and can 
play an important role in breaking the known fly-ag degeneracy from observables snch 
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as cluster number counts im. 

The evolution of the BAO peak shape is dominated by the nonlinear process of 
random displacement of galaxies with respect to their initial locations. We have shown 
that it is reliably described by the propagator part of the nonlinear PS (the hrst term in 
([^), whereas the contribution from the mode-coupling part basically disappears when 
taking ratios of the CF computed at different redshifts, or for different values of the 
neutrino masses. The propagator can be computed both in Lagrangian PT or in Eulerian 
PT after resummation (see, for instance, [211 EDI EHl EDI Ea EU), the leading result in 
both cases being the Zel’dovich propagator. The difference between the different schemes 
resides basically in the computation of the subleading corrections to the propagator and 
of the mode-coupling part, which, as we have seen, has moderate impact on the BAO 
peak shape. Therefore, any of the semi-analytic computation schemes proposed in 
the literature, based on resummations of the PT expansion in the Eulerian or in the 
Lagrangian frameworks, should give comparable performances as long as they give a 
mode-coupling independent part of the PS of the Zel’dovich form. 

We have shown that the matter CF computed from the simplest formula for the PS, 
eq. (|^, agrees with N-body simulations in real space. The effect of the same physical 
processes can also be taken into account in redshift space, leading to an agreement with 
simulations of the same quality as for the real space. 

In both cases, the quantity responsible for the BAO peak blurring by nonlinear 
effects is the velocity dispersion computed in linear theory according to eq. (0. In 
the massive neutrino cases we use the linear PS for the total matter density perturbation, 
consistently with the fact that this is the source of the gravitational potential in Poisson 
equation. 

In the case of halos our results are consistent with no velocity bias (but this topic 
should be reconsidered with larger simulations) and a nonlocal density bias of the form 
b\o + boi kf, as proposed in [371 EH l39] . 

The analysis performed in this paper for the neutrino masses can be extended to 
other cosmological parameters. In this regard, analytic formulae such as (15) can be 
extremely useful for a fast assessment of parameter forecasting for future surveys. 
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In this appendix we derive the general expression Q for the nonlinear PS. We express 
the density and velocity divergence held, Oik,z) = ik ■ v(k, 2 ;), using the compact 
notation [12] </9a(k, 2 ;), where ipiik^z) = S{k,z) and 9 ? 2 (k, 2 ;) = —6'(k, 

We assume that the held (pa(k, z) at a redshift z is a deterministic functional of 
the initial density and velocity helcjH dehned at some initial redshift Zin, and 

therefore it can be expressed as a functional expansion 


n=\ 


dfqi ■ ■ ■ d^q„ 


5^^a[ip^^]{k,z) 


TZMi)---TZS^n). 


(A.l) 


where momentum conservation implies that the coefficients of the expansion are oc 

^^(k-Er=iq*)- 

The propagator is dehned as the ensemble average of the hrst derivative of 
93 a[ 93 *”] (k, 2 ;) with respect to 93 *”, 


(<(qi)---<(qn)> . (A.2) 


1 r 

n=0 ■ 


^^+Va[93^^](k,2:) 

h^pZ{af}6^pZMi) ■ ■ ■ hq^Z('A^ 


where we have used translation and rotation invariance of the ensemble averages, that 
is (93a"(qi) • • • 95c”(qn)) oc {‘2.T^fhD{YJi=i q*) ^^^(k, z) = Gac{k, z). We note that we 

have not assumed that the initial conditions are gaussian. 

The PS is given by 


(27r)^^£,(k + k')Fab(k, z) = (93a(k, ^)93fe(k', z)). 


(A.3) 


Expanding both functions as in 
of the form 


(|Alll we get a sum of terms containing held averages 


«(qi)---<(qnVS(pi) 


■ Tdm 


(Pm)) , 


(A.4) 


where the hrst n terms come from the expansion of 93 a(k, 2 ;) and the last m ones from 
that of 93fe(k',z), and the sum over n and m is taken. We can isolate all the terms in 
which only one of the 9?c"(qi) is averaged with one of the y^Z-{pj), and we will collectively 
indicate as “mode coupling” terms all the contributions to the above averages which are 
not of this form. Eq. (A.4) can then be rewritten as 


nm (92ci(qi)<dd"(Pi))(‘dS(q2) ■ • Wc’I(qn))(<dd 2 (P 2 ) • • Wdl(Pm)) +mode coupling 
= nm {27if6D{qi + Pi)Pcdi(li){TZ (^'^)''' TZ(^n)}{TZ(P 2 ) ■ ■ ■ TZ(Pm)} 

+ mode coupling, (A.5) 


§ We assume here that vorticity and higher moments of the distribution function (i.e. velocity 
dispersion....) give negligible contribution to the scales we are interested in. However, the derivation 
can be easily extended to include initial conditions for these quantities as well, and the same form 0 
would be obtained in this case, after a proper redefinition of the propagator. 
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where (27r)^(5_D(q + p)Pff;{q) = As we see in the “mode-coupling” 

contributions the external momentum k and k' are split in two or more momenta of 
the initial fields The n m factor in front of the above expression comes from the 
symmetric integration and summation in (A.l). 


Using (A.5) and the expression for the propagator (A.2) we can write the nonlinear 
PS as 




(27r)^5£,(k -h \i!)Pah{k, z) = 

/'dy<iy(2x)%(q+p)p5(,) ^ 

= (2,r)='fo(k + k') (G„(k, z)Gm(J;, z)Pi2(k) + P^'^(k, z)) . 


+ mode coupling, 
(A.6) 


If we assume growing mode initial conditions then = </? 2 ^(q), and therefore 

Pli{q) = Pi 2 {(l) = ^in)- Moreover, if we are interested in the density 

PS we have to take a = b = 1. This gives exactly the expression in Q with 

G{k, z) = Gu{k, z) + Gu{k, z ), Pucik, z) = P,^^{k, z). (A.7) 


Appendix B. ZePdovich approximation 


In this appendix we discuss the relation between formula Q and the Zel’dovich 
approximation [52l [19], [53]. At a given time, which we omit in the following equations, 
particles have moved from the the initial (Lagrangian) positions q by the dispacement 
field ^'(q) = x — q, where x is the position of the particle which was in q initially. The 
matter density perturbation, 5(k), can then be expressed as 


5(k) = / (fx6(x.)e 


—zkx 


= J - 1 ) , 


(B.l) 


by means of the mass conservation relation (1 -|- ^(x)) d^x = dfq. Therefore, the PS is 
given by 


P{k) = / - 1) , 


(B.2) 


where r = q — q' and = ^(q) ~^(q0- The Zel’dovich approximation now amounts 
to using the linear expression for the displacement field, ^'(q) = ^ / 

Assuming (5^®®®(k) gaussian, the average of the exponent can be evaluated as (e“®*^'^^) = 
e-i/ 2 {(k-A«') 2 )^ and we can write the PS as 


where 



/(k,r) 


dfp (k ■ p)^ 
(27r)3 


cos(p ■ r)P^®®®(p). 


(B.3) 


(B.4) 


Eq. (B.3) provides the relation between the linear and nonlinear PS in the Zel’dovich 
approximation. It can be evaluated numerically and compared to simulations, see for 
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Figure Bl. CF for matter in real space, for massless neutrino cosmologies and at 
redshift z = 0. The green solid lines are obtained from (based on the PS eq. Q, 
that is from the Zel’dovich approximation in which the mode-coupling contribution is 
neglected, see eq. (B.7)), while the red dashed lines are obtained from the Zel’dovich 
approximation (based on the PS (B.3)). 


instance [53]. Expanding the second exponential in parenthesis in (B.3) we can also cast 
the PS in the following form 

cPpi cPpn 


= e 




n 

n=l 


(27r)^ (27r)" 


(27r)^(5 



x[F„(pi,---,p„)fP^*"(pi)---P^*"(p„), 
where the kernel functions are given by 

1 k ■ Pi k ■ p„ 




n\ p1 


Pt, 


Isolating the hrst term in the sum ( |B.5 ), it can be rewritten as 

pzeid^k) = 


MC (^) ) 


(B.5) 


(B.6) 


(B.7) 


showing explicitly that our eq. (|^ amounts to neglecting the mode-coupling contribution 
in the Zel’dovich approximation. 

The comparison between the CF’s obtained from eq. (|^ and from the full Zel’dovich 
approximation, eq. (B.3), is shown in Figures [ bI| and B2, from which one sees clearly that 
the mode-coupling part is completely negligible in the BAO peak region. Moreover, we 
also see that the two functions give essentially an identical result for ratios of correlation 
functions (both between different redshifts within one cosmology, and between different 
cosmologies with different neutrino masses). 


Appendix C. Redshift space 

We want to describe the effect of bulk flows on the CF or on the PS, in redshift space. 
In real space, the position of a particle at time r is given by 

x(r) = x(0) -b f (ir'v(r') . (C.l) 

Jo 
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Figure B2. Ratios of matter CFs in real space. On the left panel we show ratios 
of CFs at different redshifts, for massless neutrinos. On the right panel we show the 
z = 0 OF for massive neutrinos (the label on the curve giving ffmi, in eV), divided 
by the corresponding CF for massless neutrinos. The green solid lines are obtained 
from (based on the PS eq. Q), while the red dashed lines are obtained from the 
Zeldovich approximation (based on the PS (B.3)). 


In redshift space it is shifted to 

v^(t) 

s(r) = x(r) + z , (C.2) 

where TL is the Hubble parameter computed using conformal time, we work in the distant 
observer approximation and we set the line of sight along the z-axis. Now we split the 
particle velocity into a short and a long distance component, 

v(r) = ^sh{,r) + v;(r), (C.3) 


so that 


where 


s(r)=y(r) + ^z + ^ + ^fh),. 


n 


nf n 


y(r) = x(0) + / dr'v^fe(r'), 

Jo 

and we have assumed that v/ evolves according to linear PT. 
Then, imposing mass conservation 

(1 + 5sis))£s = (1 + 5{y))d^y , 


(C.4) 

(C.5) 


(C.6) 


and Fourier transforming, we have the density contrast in redshift space 

^^(k) + 5,(k) = J (1 + 5(x)). (C.7) 

The statistical average over long modes only gives 

^^(k) + {6s{k))i = J (1 + 

f flc ^ /(2 + /)) I'll ‘ ^ ' ...... 

= / d-^xe **"^6 ^ e (1 + (5(x)), (C.8) 


where /i = k ■ z. 
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Figure Dl. Two-point correlation function of the matter density field measured 
using 10 N-body simulations with 256^ (blue), 384^ (green) and 512^ (red) particles. 
The error bars represent the scatter around the mean value for the 10 realizations 
with 256^ particles divided by the square root of 10. The bottom panel displays the 
different correlation functions normalized to the fiducial one (256^). 


We now expand the third exponential in (C.8) to linear order in Kh and then we 
get the correlator 

{{i.(k)),(i,(-k)),).i = , (C.9) 

where 




Kaiser 

sh 


(F/i) = (i+/tvmr(fc), 


(C.IO) 


with Pfff{k) the linear PS for the short modes. 


Appendix D. Resolution tests 

Here we study the dependence of our results with respect to mass and force resolution. 
In order to check whether our results are converged against mass and force resolution, we 
have run two different sets of simulations for the massless neutrino cosmological model; 
10 simulations with 384^ CDM particles and 10 simulations with 512^ CDM particles. 
All simulations have a box size equal to 1000 Mpc/h and the softening length is set to 
1/40 of the mean linear inter-particle separation. 

We have then measured the matter two-point correlation function in the different 
simulation sets and we show the results in Fig. 


Dl 


In that hgure we also show the 
results obtained by using 10 realizations of our fiducial suite containing 256^ particles 
with the error bars representing the scatter around the mean value divided by vTo. 

We hnd that by increasing the simulation mass and force resolution the differences 
in the matter correlation function, in the range of scales investigated in this paper, are 
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